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INTRODUCTION 


Exponentially  weighted  Galerkin-f inite  element , ^ ^ ^  collocation,^  and 
exponentially  fitted  finite  dif f erence1 » 5 » 7  schemes  have  become  popular  and 
effective  numerical  methods  for  solving  convection  dominated  convection- 
diffusion  problems.  They  avoid  the  spurious  mesh  oscillations  found  in 
centered  schemes  at  high  values  of  the  cell  Reynolds  or  Peclet  numbers  and 
reduce  the  effects  of  numerical  diffusion  found  in  upwind  finite  difference 
schemes. 

The  exponential  schemes  all  require  evaluating  the  function 

E,  =  f(z)  :=  coth  z  -  1/z 

in  order  to  obtain  their  optimal  accuracy.  For  example,  the  exponentially 
fitted  Galerkin-f inite  element  method  for  the  two-point  boundary  value  problem 

s  .  c(x)  --  =  0  ,  0  <  x  <  1  ,  u(0)  =  A,  u(l)  =  B  (2) 

dyr  dx 

on  a  uniform  grid  of  spacing  h  =  1/N  is  given  by,  (cf ,  e.g.  ,  Hughes  ) 

(e/h2  ) ( U-l—  !  -  2Ui  +  Ui+1)  -  ( l/2h)  [  ( 1+Ci)c(xi)  (UfUi- L) 

+  (l-5i)c(xi+i)(Ui+1-Ui)]  =  0  ,  i  =  1,2,... ,N-1  (3a) 

UQ  =  A  ,  UN  =  B  (3b, c) 

Here  denotes  the  numerical  approximation  of  u(ih),  i  =  0,1,..., N,  x^  is 
some  point  on  (xi-i,  xi) ,  e.g.,  the  center  of  the  subinterval,  and  can  be 


*Ref erences  are  listed  at  the  end  of  this  report. 
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interpreted  as  a  function  evaluation  point  in  a  one-point  quadrature  rule  (cf. 


Hughes^).  The  choice 


5i  -  f(Pi/2)  (4) 

where  Pi  is  the  cell  Reynolds  or  Peclet  number 

Pi  =  c(xi)h/ e  (5) 

is  known  to  give  the  exact  solution  of  Eq.  (2)  for  all  pi  when  c  is  a  constant 
(cf;  Christie  et  al.^). 

The  function  f(z)  is  relatively  expensive  to  evaluate  because  of  the 
exponential  functions  and  is  usually  replaced  by  the  "doubly  asymptotic" 
approximat ion 

z/3  ,  I z |  <3 

5  =  f a( z)  :=  (6) 

sgn(z),  | z |  >  3 


The  function  f^(z)  provides  an  0(z3)  approximation  to  f(z)  when  |z|  <<  1  and 
an  0(l/z)  approximation  when  |z|  »  1.  Furthermore,  when  Eq,  (6)  is  used  in 
Eq.  (3)  and  c(x)  is  smooth,  provides  an  0(p4)  approximation  when 

p :=  max  I p ^ I  (7) 

l<i<N 

is  small  and  an  0( 1/p )  approximation  when  p  is  large. 

Thus,  f^(z)  provides  a  good  approximation  of  f(z)  when  z  is  either  small 
or  large,  but  has  large  errors  when  z  -  0(1)  (cf.  Figure  1).  The  largest 
difference  between  f(z)  and  f^(z)  is  0.328  and  it  occurs  at  z  =  3.  This 
corresponds  to  a  value  of  p  *  6  and  since  cell  Reynolds  numbers  in  this 
vicinity  are  reasonably  common  in  computation  it  behooves  us  to  find  a  better 
approximation  for  f(z)  than  f^(z)  when  z  »  0(1). 


2 


In  this  note,  we  consider  rational  function  approximations  having  the 


form 


£  =  f r(  z  )  :  = 


z  ( 1  +  a|z| ) 

3  +  3 1 z |  +  az2 


(8' 


for  appropriate  choices  of  a  and  3*  This  approximation  will  be  considered 
successful  if  it  provides  better  accuracy  than  f^Cz)  and  is  still  less 
expensive  to  evaluate  than  f(z). 

Like  f^(z),  we  see  that  fR(z)  correctly  approximates  f(z)  as  z  +  0  and  as 
|z|  +  00  for  all  values  of  B  and  all  a  t  0.  The  maximum  difference  between 
f(z)  and  f  r( z )  is  about  0.0115  for  the  nearly  optimal  values  of  a  =  0.6  and  3 
=  1.38  (cf.  Figure  1).  Furthermore,  when  f,  f^,  and  fR  were  evaluated  for 
1000  values  of  z  e  [0,100]  we  found  that  fR  took  35  percent  less  time  to 
evaljate  than  f  while  f^  took  49  percent  less  time  than  f.  The  approximation 
fR  also  provided  greater  accuracy  than  f^  for  the  computed  solution  U^,  i  = 
0,1,..., N,  of  two  model  problems.  The  savings  in  time  and  improvement  in 
accuracy  are  significant  and  may  be  especially  important  in  multi-dimensional 
problems.  As  previously  noted,  the  greatest  gains  occur  when  =  0(1)  and 
c(x)  is  smooth. 


RATIONAL  FUNCTION  APPROXIMATION 


We  will  want  to  restrict  a  and  3  in  Eq.  (8)  so  that  our  approximation 
fft(z)  satisfies  the  following  three  conditions: 

i.  fR(z)  should  be  a  good  approximation  of  f(z)  when  z  is  small  and 
large,  although  not  necessarily  as  good  as  f^(z). 
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Figure  1.  Optimal „  doubly  asymptotic,  rational,  and  critical  choices  of 
K  =  f(z)  given  by  equations  (1),  (6),  (8),  and  £  =  1  -  1/z, 
respectively. 
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ii.  The  solution  of  Eq.  (3)  should  be  oscillation  free.  Christie  et 
al.4  have  shown  that  this  will  be  the  case  when 

z  >  0  ,  5  >  1  -  1/z 

z  <  0  ,  5  <  -(1  -  1/z)  (9) 

z  =  0  ,  all  5 

iii.  Since  f(z)  is  a  monotonically  increasing  function  of  z  we  ask  that 

dfR(z) 

— - -  >  o  ,  all  z  (10) 

dz 


S_nce  both  f(z)  and  fR(z)  are  odd  functions  of  z  it  suffices  to  enforce 
these  conditions  for  z  >  0.  We  shall  see  that  enough  flexibility  remains  for 
us  to  select  a  and  £  to  improve  accuracy  when  z  =  0(1). 

The  functions  f(z)  and  fR(z)  have  the  following  asymptotic  behavior  for 
small  and  large  values  of  z: 


z  z3  2z5 

0(z7)  ,  z  «  1 

3  45  945 

f(z)  = 

1  -  -  +  0(e~2z)  }  z  »  1 

z 


(11a) 


(lib) 


z  1  11 

-  [1  +  (a  -  -B)z  +  -  (-32-a3~a)z2  +  0(z3)],  z  «  1  (12a) 

3-1  32-3-3cx  1 

fR(z)  =  1 - +  —5-5—  +  0(-~)  ,  z  »  1  ,  a  *  0  (12b) 

az  a4z4  z5 

3  9  1 

t1  ~  ~~  +  ~2~Z  +  °(~3^  »  z  »  1,  a  =  0  ,  3*0  (12c) 

3z  3  z4  z° 


1 

3 
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Equation  (12c)  only  gives  the  correct  limiting  value  of  f(z)  as  z  -►  00  when  3  = 
1  and  this  value  of  3  does  not  satisfy  Eq.  (9);  hence,  we  will  no  longer 
consider  approximations  with  a  =  0. 

Let 


e(z)  :=  fR(z)  -  f ( z ) 

denote  the  pointwise  error  and  use  Eqs.  (11)  and  (12)  to  obtain 


(13) 


1  z2  i  l  zi 

(a  -  -3)"  +  (-32-cx3-a+-)~  +  0(z4) 
3  3  3  5  9 


z  «  1 


e(z)  = 


(14) 


3-1  1  32-3-3a  1 

(1 - )-  +  —3-5-  +  0(~|) 

a  z  orz^  zJ 


z  »  1 


We  see  that  the  rate  of  convergence  as  z  -*•  0  can  be  improved  from  0(z2)  to 
O(z^)  by  selecting 

3  =  3a  (15) 

while  the  rate  of  convergence  as  z  -*•  00  can  be  improved  from  0(l/z)  to  0(l/z2) 
by  selecting 

3  =  1  +  a  (16) 

Both  Eqs.  (15)  and  (16)  can  be  satisfied  simultaneously  by  selecting  a  =  1/2 
and  3  =  3/2. 

Before  deciding  on  either  or  both  of  Eqs.  (15)  and  (16)  we  still  must 
find  bounds  on  a  and  3  so  that  Eqs.  (9)  and  (10)  are  satisfied.  It  is 
slightly  simpler  to  consider  Eq.  (10)  first;  thus,  we  differentiate  Eq.  (8)  to 
obtain 

dfR  3  +  6az  +  az2(3~l) 

—  =  — - - -  >  z  >  0  <17) 

dz  (3+3z+az-)2 

For  dfft/dz  >  0,  the  polynomial 

p(z)  =  3  +  6az  +  az2(3~l)  (18) 


6 


shoulc  not  have  any  positive  roots.  It  will  have  two  negative  roots  if  a  >  0 
and  8  >  1  and  two  complex  roots  if  a  >  0  and  8  >  I4*3a  or  a  <  0  and  8  <  14-3a. 
For  the  reasons  of  accuracy  expressed  by  Eqs.  (15)  and  (16)  we  would  like  to 
be  as  close  to  a  =  1/2  and  8  =  3/2  as  possible.  Hence,  we  will  not  consider 
the  region  where  a  <  0  and  confine  our  attention  to  choices  satisfying  a  >  0 
and  8  >  1* 

•  Finally  using  Eq.  (8),  condition  (9)  will  be  satisfied  if 


1  3  4*  (8“3)z  4*  (l4*a-8)z2 

fR(z)  -  (1  -  -)  - - 57 - 5; -  >  0 

z  z2  (34*8z4*az2  ) 


Since  a  and  8  are  positive  we  want  the  polynomial 

p  ( z )  =  3  +  ( f?-3)z  +  (l+a-e)z2 
to  have  no  positive  roots.  p(z)  will  have  two  negative 
and  two  complex  roots  if 


,  z  >  0  (19) 

(20) 

roots  if  3  <  8  <  I4*a 


8  <  -3  4-  2/3(l+a) 


(21) 


The  values  of  a  and  8  that  satisfy  both  Eqs.  (9)  and  (10)  are 


-  <  a  <  3  ,  1  <  8  <  -3  4-  2/30+a) 

3 

a  >  3  ,  1  <  8  <  1+a  (22) 

This  region  is  shown  shaded  in  Figure  2.  Note  that  the  point  a  =  1/2,  8  = 

3/2,  which  improves  accuracy  for  small  and  large  values  of  z,  fails  to  satisfy 
condition  (10).  However,  Figure  2  suggests  that  an  effective  alternative 
might  be  to  pick  the  point  on  the  curve  8  =  '3  4-  2/3(  14-a)  that  is  closest  to  a 
=  1/2,  8  =  3/2.  This  point  is  a  =  0.5931,  8  =  1.3723  and  a  search  shows  that 
it  is  near  the  point  which  minimizes  the  maximum  value  of  |e(z)|,  for  all  z. 
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Figure  2.  Region  of  acceptable  values  of  a  and  3  (shaded)  that 
satisfies  equations  (9)  and  (10) . 
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In  the  examples  of  the  next  section  we  used  a  =  0.6  and  B  =  1.38  for  which 
the  maximum  value  of  |e(z)|  is  0.0115. 

NUMERICAL  RESULTS 

The  approximations  f^(z)  and  fR(z)  and  the  exact  function  f(z)  given  by 
Eqs .  (6),  (8),  and  (1),  respectively,  were  coded  in  FORTRAN  and  used  in  the 
series  of  numerical  experiments  described  below. 

Equation  (1)  cannot  be  used  for  f(z)  when  z  is  near  zero  and  it  was 
replaced  by  the  expansion  (11a)  for  |z|  <  0.01.  The  FORTRAN  Library  function 
for  tanh  (z)  was  used  in  Eq.  (1)  when  z  was  not  small.  All  calculations  were 
performed  in  double  precision  arithmetic  on  either  an  IBM  3033  at  the 
Rensselaer  Polytechnic  Institute  or  an  IBM  4341  at  the  Benet  Weapons 
Laboratory . 

In  our  first  experiment,  we  evaluated  f(z),  f^(z),  and  fft(z)  for  1000 
values  of  z  £  [0,100]  and  timed  the  results.  The  normalized  times  recorded  in 
Table  1  were  averaged  over  several  runs  and  include  only  the  times  to  evaluate 
the  functions  and  neither  input/output  nor  supervisor  state  times.  Variations 
in  times  from  run  to  run  were  less  than  two  percent.  The  results  indicate 
that  f r  took  35  percent  and  f/^  took  49  percent  less  time  to  evaluate  than  f. 
The  relative  timing  figures  can  be  expected  to  vary  significantly  from 
computer  to  computer  and  even  from  compiler  to  compiler;  however,  the 
differences  between  the  times  to  evaluate  fR  and  f  are  large  enough  so  that 
savings  should  be  achieved  in  most  environments. 
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TABLE  1.  RELATIVE  COMPUTER  TIMES  TO  EVALUATE  THE  APPROXIMATIONS  fA(z) 


AND  fR(z)  AND  THE  EXACT  FUNCTION  f(z)  FOR  1000  VALUES  OF  z 
ON  0  <  z  <  100. 


Method 


Time 


1 

Doubly  Asymptotic 

Approximation,  Eq.  (6) 

1 

0.508 

i 

Rational  Approximation, 

Eq.  (8)  1 

0.654 

1 

Exact,  Eqs.  (1),  (11a) 

1 

1.000 

TABLE  2.  MAXIMUM  ERRORS  AT  STEADY  STATE  FOR  EXAMPLE  1. 


1 

1 

i 

1  Doubly 

1 

Lax 

1 

1 

T 

1 

P 

1  x 

1 

I  Asymptotic 

1 

Wendrof f 

i 

|  Rational 

i 

1 

1 

Optimal 

6 

T~ 

I  0.75 

1 

1 

|  0. 204x10“ 1 

1 

i 

1  0. 159xl0“2 

1 

1  0.866x10“ 3 

1 

1 

1 

l 

0.250x10“ 1 3 

500 

I  0.95 

i 

|  0. 200x10“ 2 

1 

1 

I  0.235x10“ 1 

1 

I  0.705x10“ 3 

1 

1 

1 

1 

0.208x10“ 1 1 
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In  order  to  study  variations  in  the  computed  solutions  when  the  different 
approximations  of  f(z)  are  used  we  consider  two  boundary  value  problems  having 
the  form 


d2u  dF(u) 

e  - =0  ,  0  <  x  <  1  ,  u(0)  =  A  ,  u(  1 )  =  0  (23) 

dx2  dx 


Since  our  motivation  for  performing  the  work  described  in  this  note  is  to 
study  exponentially  weighted  finite  element  schemes  for  transient  problems  we 
do  not  use  the  numerical  method  (3)  but  rather,  we  follow  Osher^  and  consider 
■the  solution  of  (23)  as  the  steady  state  limit  of  the  following  initial- 
boundary  value  problem 

3u  3F(u)  3^u 

~  + - =  e  — -  ,  0  <  x  <  1  ,  t  >  0  (24a) 

3t  3x  3xz 


A,  x  =  0 

u(0,t)  *■  A  ,  u(l,t)  =  0  ,  u(x,0)  =  (24b, c,d) 

0,  0  <  x  <  1 

Equations  (24)  are  approximated  by  the  following  explicit  difference  scheme 


n+1  n  n  n  n  n 

Ui  =  Ui  -  X  [(l+5i)(Fi-Fi.1)  +  (l-?i)(Fi+1-Fi) 


n  n  n 

+  (eA/h)(Ui-1-2Ui+Ui+1)  ,  i  =  1,2,...,N-1  (25a) 


n  n  o  x  A  w 

U0  »  A  ,  UN  =  0  ,  Ui  *  (25b, c,d) 

0,  i  =  1,2,.. • ,N 

n 

where  Ui  denotes  the  numerical  approximation  of  u(ih,nAt),  At  is  the  time 
step,  and 

X  =  At/h  (26) 
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The  cell  Reynolds  number  is  still  given  by  Eq.  (5)  with  c  defined  as 

dF(u(x,t) ) 

c(x,t)  - - - - 

du 


(27a) 


and 

n  n 

xi  >  ui-l  +  ui  <  0 

xi  =  (xi-i+Xi)/2  ,  u"_i  +  uj  =  0  (27b) 

xi-l  »  (ui-l+ui)  >  0 


For  reasons  discussed  in  Flaherty  and  Hathon4  these  choices  of  xi  might  be 
better  than  always  using  the  center  of  the  subinterval  (xi-^x^). 

The  explicit  difference  scheme  will  be  stable  to  linear  perturbations 

provided  that  (cf.  Osher8  or  Flaherty3) 

A[c(x^,nAt)  +  2e/h]  <1  ,  i  =  1,2,...,N,  n  >  0  (2 


We  use 


E  :=  max|u(ih,nAt)  -  U-£  | 
l<i<N 


(29) 


th  n  chosen  large  enough  so  that  steady  state  has  been  reached,  to  measure 


wi 


errors • 


Example  1:  Consider  the  constant  coefficient  problem  for  Eq.  (23)  with 
F(u)  =u,  c  =  1,  A  =  1 ,  which  has  the  exact  solution 


1  -  e 


-(l-x)/e 


u(x)  = 


(30) 


-  p-l/e 


1  -  e 
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The  results  of  calculations  when  Si  in  Eq.  (25)  was  evaluated  by  the  doubj.. 
asymptotic  approximation  (6),  the  rational  approximation  (8),  the  Lax-Wenc -o 
scheme  (Si  =  X),  and  the  optimal  scheme  (1)  are  presented  in  Table  2  for  N  - 

20,  p  -  6,  X  -  0.75,  and  N  -  20,  p  -  500,  X  -  0.95.  The  optimal  scheme  (1)  is 

exact  for  this  example.  The  samll  errors  reported  in  Table  2  are  due  t ^ 
combined  effects  of  roundoff  and  our  failure  to  reach  the  steady  state  limit. 
As  expected  the  rational  approximation  improves  upon  the  results  of  the  doubly 

asymptotic  approximation  and  the  improvement  is  greatest  for  p  -  6.  The 

Lax-Wendroff  solution  oscillates  when  x  is  near  unity. 

Example  2.  We  consider  the  nonlinear  Burgers'  equation,  F(u)  -  u2/2 
c  =*  u,  A  ■  tanh  l/2e,  for  which  the  exact  solution  of  Eq.  (23)  is 

u(x)  -  tanh(l-x)/2e  (31) 

Results  comparing  the  doubly  asymptotic,  rational,  and  optimal  choices  of 
are  shown  in  Table  3  for  N  ”  20,  p  ”  h/e  m  6,  X  =  0.75,  and  N  •*  ^.0,  p  *  500,  X 
-  0.95,  and  in  Table  4  for  e  -  1/128,  p  -  1,2,4,8,16.  In  Table  5  we  show 
results  for  the  error  |Ut  -  u(ih,nAt)|  at  steady  state  and  x  =  ih  -  0.875, 
0.9375  for  e  =  1/128  and  p  =  1,2,4,8,16.  The  solution  obtained  using  the 
doubly  asymptotic  approximation  had  mesh  oscillations  for  P  ■*  8  and  X  ”  0.96 
so  this  calculation  was  rerun  with  X  -  0.48. 

Although  accuracy  is  not  as  good  for  this  nonlinear  example,  the  results 
generally  parallel  our  findings  for  the  linear  problem.  Table  5  shows  that 
the  rational  and  optimal  choices  of  are  better  than  the  doubly  asymptotic 
choice  at  reducing  the  effects  of  numerical  diffusion  for  cell  Reynolds 
numbers  in  the  range  of  2  to  16. 
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TABLE  3.  MAXIMUM  ERRORS  AT  STEADY  STATE  FOR  EXAMPLE  2. 


'T 


6 

500 


0.75 

0.95 


Doubly 

Asymptotic 


0.124 


0.200x10 


-2 


Rational 


0.761x10' 1 
0.135xl0'2 


TABLE  4.  MAXIMUM  ERRORS  AT  STEADY  STATE  EUR  EXAMPLE 


AN  *  DENOTES  THAT  X  =  0.48 

FOR  THIS  CASK. 

1 

1 

1 

1 

1 

Doubly 

7 

1 

1 

1 

p  1 

1 

x  1 

1 

Asymptotic 

Rational 

j  _ 

1 

1 

16  | 

""  r 

0.96  | 

1 

0.605x10' 1 

1 

|  0. 330x10' 1 

1 

1 

1 

1 

1 

8  1 
i 

0.96  | 

I 

0.117* 

1 

|  0. 599x10' 1 

1 

1 

1 

1 

1 

4  I 

i 

1 

0.96  | 

1 

0.110 

1 

|  0.940x10' 1 

1 

1 

1 

2  1 

1 

0.6912  | 

1 

0.627x10' 1 

1 

|  0.627x10' 1 

1 

1 

1 

1 

1  | 

1 

0.4608  | 

1 

0. 156x10' 1 

1 

I  0. 157x10' 1 

1 

Opt imal 

i 

0.7  66x10' 1  | 

I 

O.lOOxlO'2  I 

_ _ _ L 

2  WITH  e  =  1/128. 

t - r 

i  i 

|  Optimal  | 

I _ | 

|  0. 308x10' 1  I 

I  I 

|  0. 601x10' 1  | 

I  I 

|  0.932x10' 1  I 

I  I 

|  0.618x10' 1  | 

I  I 

|  0.1 56x10' 1  | 
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TABLE  5.  ERRORS  AT  =  0.875  AND  0.9375  AT  STEADY  STATE  FOR  SAMPLE  2  WITH 
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CONCUIS  [ON 


We  have  shown  that  the  rational  function  approximation  (8)  is  an 
alternative  to  the  doubly  asymptotic  approximation  (6)  of  f(z)  that  offers 
greater  accuracy  for  about  a  30  percent  increase  in  cost.  The  approximation 
is  most  useful  for  cell  Reynolds  numbers  in  the  range  of  one  to  ten. 
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